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Abstract 

We study chaotic systems with multiple time delays that range over several orders of magnitude. 
We show that the spectrum of Lyapunov exponents (LE) in such systems possesses a hierarchical 
structure, with different parts scaling with the different delays. This leads to different types of 
chaos, depending on the scaling of the maximal LE. Our results are relevant, in particular, for the 
synchronization properties of hierarchical networks (networks of networks) where the nodes of sub- 
networks are coupled with shorter delays and couplings between different subnetworks are realized 
with longer delay times. Units within a subnetwork can synchronize if the maximal exponent scales 
with the shorter delay, long range synchronization between different subnetworks is only possible 
if the maximal exponent scales with the long delay. The results are illustrated analytically for 
Bernoulli maps and numerically for tent maps. 

PACS numbers: 05.45.Xt, 89.75.Hc, 02.30.Ks 
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Networks of nonlinear units with time-delayed interactions play an important role in vari- 
ous systems, such as coupled semiconductor lasers, predator/prey systems, traffic dynamics, 
communication networks, genetic circuits, or the brain [l-6]. Delay times may induce high- 
dimensional chaotic dynamics [l-9] as for example in semiconductor lasers with delayed 
feedback Q. A p artlCularly interest phe_ in this eontext is the zero lag synchro- 



11 



151 ] . Chaos synchronization 



nization of chaotic units, despite the long interaction delays 
finds applications in encrypted communication fli]. \l7\. 

Chaos in the network is quantified by the spectrum of Lyapunov exponents (LEs), which 
measures the sensitivity to initial conditions. For a system with one long delay, chaos can be 
characterized by the scaling of the maximal LE with increasing delay time r. In the region 
of strong chaos, it approaches a positive constant value whereas for weak chaos the maximal 
LE decreases as 1/r. These scaling properties have consequences for chaos synchronization: 
Networks with strong chaos cannot synchronize completely for long delays, whereas for weak 
chaos synchronization is possible depending on the value of the maximal LE and the topology 
of the network Q, 19| . 

Strong and weak chaos have been demonstrated in networks with a single time delay, up 
to now. However, realistic systems may have different transmission delays for the coupling 
signals. In a network with a distribution of delays complex behavior is expected to be 
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J-I22J . However, if the network has different delay times with special integer 
ratios, one finds resonances which can either stabilize or rule out chaos synchronization 23], 
depending on the ratio. 

In this Letter we study networks of nonlinear units coupled by multiple delay times which 
differ by several orders of magnitude. A typical example is a network of networks, with a 
connection delay T\ between the nodes within a sub-network and a much longer connection 
delay T2 between the different sub-networks. We explain the scaling of the full spectrum 
of LEs with increasing delay times and extend the concepts of strong and weak chaos to 
multiple delay systems. Finally, we relate the synchronization properties of a hierarchical 
network to the scaling behavior of the LEs. 

Spectrum of Lyapunov exponents of a Bernoulli map. Strong and weak chaos have been 
found both for time-continuous and discrete sytems with delay js, (2^. Since the main 
results are valid in both cases, we perform our calculations for iterated maps. For networks 
of Bernoulli maps, even analytic results can be derived. We start with a single chaotic map 
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with N different feedback delays t±, . . . , r^r, described by 

N 

x t +i = (1 - e)f(x t ) + e^\ fc /(x t _ T J , (1) 

k=i 

with xel. We consider here delays with different orders of magnitudes 1 <C T\ <C . . . <C Tn- 
The spectrum of LEs A = {Ai, . . . , A^}, which describes the evolution of a perturbation 5x t 
along a chaotic trajectory s t , is calculated using the linearized equation 

N 

Sxt+i = (1 - e)f(s t )Sx t + ^kf'{s t ^ Tk )5x t - Tk . (2) 

fc=i 

The coefficients f'(st) in general depend on time, since the trajectory s t is time dependent. 
However, for a Bernoulli map, given by 

f(x) = ax mod 1 (3) 

with \a\ > 1, we have constant coefficients f'(x) = a. The linearized equation ([2]) reduces 
then to a polynomial equation for the characteristic multipliers z, which characterize the 
growth of a perturbation 5x t = SxqZ 1 : 

N 

z = (l-e)f' + ef'J2^z- Tk . (4) 

k=l 

The Lyapunov exponents A e A are given by A = In \z\. The calculation of the LEs can thus 
be performed in the same way as the stability calculation of a steady state. For a single delay 
system, it is known that such an equation can have two different types of unstable (with 
\z\ > 1) solutions 2|| : strongly unstable and weakly unstable. The strongly unstable 
root is approximated by the delay-independent term z w (1 — e)f, provided |(1 — e)f'\ > 1, 
and does not depend on the delays to the leading order. Also the multiple-delay system can 
have a strongly unstable multiplier. 

In analogy with the single delay system, we assume a scaling behavior for the next group 
of multipliers Inz = ioj + 7i/ti, with ji > 0. In leading order, by neglecting all terms of 
order l/r±, e~ 7ir2 / ri , and smaller, the characteristic polynomial 01]) reads 

e iu = (1 - e )f + eKj'e-^ 1 '^ . (5) 

This equation allows to compute a curve 71 (u) = In \ eKif'\ — In \ e tuJ — (1 — e)/'|, on which 
the roots are located. The imaginary parts of these roots differ by approximately Au w 
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FIG. 1: (Color online) The spectrum of LEs of a Bernoulli map (Eq. (T4J) (red dots) subject 
to three feedback delays. The different panels are zooms. The full black lines are the analytic 
approximations for long delays 7i(w)/ti (panel (b)), 72(w)/t2 (panel (c)) and 73(w)/t3 (panel (d)). 
Parameters are a = 3, e = 0.63, n\ = 0.3, K2 = 0.6, K3 = 0.1, t\ = 40, T2 = 500 and T3 = 6000. 

27r/ri, the number of roots thus increases linearly with 77. In the following we call the 
unstable subset of these LEs the Ti-spectrum A x . The Ti-spectrum corresponds to the 
pseudo-continuous spectrum for steady states of single-delay systems [25]. 

For the third group of unstable multipliers we assume a scaling with the second delay 
\nz = iu + 72/V2, with 72 > 0. In leading order, we find a characteristic polynomial 



/ —lUJTl 



-10JT2 +72 



(6) 



In an analogous way, the unstable roots (if they exist) are approximated by a curve 72 (00) = 
In \eK2f'\ — In \e lul — (1 — e)f — e~ %WTL eKif'\. The corresponding LEs form a T2-spectrum A 2 , 
which scales inversely with the second delay r 2 . The number of exponents in this spectrum 
scales linearly with t^. 

It is possible to calculate unstable spectra related to each of the delays 77; only the tn- 
spectrum, related to the largest delay present in the system, can have a stable part. Figure 
[1] shows the spectrum of LEs for a Bernoulli map, obtained by solving Eq. (j4j) numerically, 
and the analytical long delay approximations 7fc(w). Although the different time scales are 
not so far apart, the analytical curves 7fe(w) provide a good approximation for the LEs. We 
can clearly distinguish one strongly unstable multiplier, and the weakly unstable spectra 
with their respective delay-scaling. 

A similar hierarchy of eigenvalues for steady states of time-continuous delay differential 
equations with multiple delays can be shown applying the same arguments to the corre- 
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sponding transcendental equation. 

General case. The slope of a chaotic map in general depends on the trajectory, so that 
the linearization (Eq. (J2])) is time-dependent. The Lyapunov spectrum is then evaluated 



numerically using a Gram-Schmidt orthogonalization procedure according to Farmer [26]. 
We find that the properties which we derived for the time-independent case are preserved 
in the presence of fluctuations, so that the complete Lyapunov spectrum A is composition 
of Tfc-spectra. The r^-spectrum is obtained numerically by integrating the evolution of an 
auxiliary perturbation variable 8x*, for which the delay terms r^+i to tn are removed 

8x° t+1 = (1 - e)f( St )6x° t 
A = {A } 

5xl +1 = (1 - e)f(s t )5xl + e^f (8^)5x1^ (7) 

=> Ax = {Ai ,niax7 • • • j 

Ai,„}\ A 



The first exponent Ao is called the sub-LE. If a partial spectrum A^ contains positive expo- 
nents, they can be said to 'survive' the introduction of further time scale separated delay 
terms, because the contribution of the additional terms becomes exponentially small for 
the corresponding Lyapunov modes. We exclude these exponents from the definition of the 
succeeding spectra, so that the partial spectra A^ do not overlap, and each spectrum A& 
scales only inversely with Tfe. 

We demonstrate the composition of the LE spectrum into different partial spectra for a 
tent map, 

{- x for < x < a 

(8) 
j^- (1 — x) for a < x < 1 

subject to two different delayed feedbacks. Figure E](a,b) compares the sub-LE Ao, the 
partial spectrum Ai and the total spectrum A. We find that the maximal LE A max G A is 
well approximated by the sub-LE Ao for a long enough delay T\. Since A depends on the 
trajectory s t , it also depends indirectly on the delays Tfe. Nevertheless, this dependence is 
negligible, as has also been reported for the Lang-Kobayashi model representing a single 



delayed feedback system [18|, [l9|. The next exponents A 2 , . . . , A 7 of the full spectrum are 



approximated by the t\— spectrum Ai; the full spectrum A deviates from Ai as the latter 
becomes negative. The second LE X 2 , which coincides with Ai max , decreases with t%. 
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FIG. 2: (Color online) Panel (a) shows Ao and the first 15 exponents of the full and the partial 
Lyapunov spectra, A and Ai, for a tent map with two feedbacks t\ = 60 and T2 = 500 in the strong 
chaos regime. Panel (b) shows the sub-exponent Ao (green 'x' symbols), two maximal LE of the 
full spectrum A (red circles), and the maximal exponent of the ri-spectrum Ai )max (blue crosses) in 
the same setup for different delays t\. In panel (c) and (d) the LE A&, k = 2, . . . , 100 for different 
values of the feedback delays are shown. Panel (c) shows the scaling of the positive LEs with n for 
fixed T2 = 500. Panel (d) shows the T2-scaling of the smaller exponents for fixed t\ = 30. Other 
parameters are a = 0.4, e = 0.4, k\ = 0.8, and K2 = 0.2 



The different scaling behaviors of the LE spectrum are depicted in Fig. [2]^c,d). The partial 
spectrum Ai as a whole (and thus the upper part of the full spectrum) scales inversely with 
the delay T\\ the number of exponents in the partial spectrum however scales linearly with 
T\. These two effects are demonstrated by plotting A^ti vs. kjr\ (with k the ranking of 
the exponent): For different delays T\ the spectrum converges to a curve for all exponents 
from the Ti-spectrum. The curves diverge for smaller exponents kjr\ > 0.2, as these scale 
with the largest delay Ti. Similar, for varying T2 the spectrum converges to a curve for all 
exponents A& from A 2 when plotting \^t 2 vs. k/r2- 

Tk-chaos. Apart from a hierarchical Lyapunov spectrum as described above, the different 
time scales can also manifest in the maximal LE A max , which can be considered as the most 
important quantity describing a chaotic system. If A > we speak of strong chaos; the 
maximal LE A max ~ Ao and does not vary with any of the delays. In the weakly chaotic 
regimes Ao < holds. If Ai jmax > we speak of Ti-chaos and A max oc \jr\. If Ai jmax < the 
second delay dominates and A max oc l/r 2 . Consequently we define Tfc-chaos as the scaling of 
A max with l/7-fc. 
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FIG. 3: Evolution of a small perturbation for a tent map with two delayed feedbacks, for (a) strong 
chaos (e = 0.35, k = 0.8), (b) n-chaos (e = 0.7, « = 0.8) and (c) r 2 -chaos (e = 0.7, « = 0.2). Other 
parameters are T\ = 47, t 2 = 500 and a = 0.4 

The difference between strong, n-, T2-,. ■ ■ , 7fc-chaotic dynamics can be directly observed 
in the evolution of a small perturbation. Fig. [3] shows the difference 5xt between trajectories 
of two identical chaotic tent maps, initialized identically except for a point-like perturbation 
at t — 0. The instantaneous evolution of this perturbation is governed by Ao- In the strongly 
chaotic regime, it thus increases exponentially, as shown in Fig. [3(a). In the weakly chaotic 
regimes the perturbation decays first, to reappear and decay at t ~ r±. The perturbation 
evolves over the consecutive Ti-intervals according to Ai. max . In the Ti-chaotic regime a 
perturbation thus spreads on the time scale of t\ (illustrated in Fig. [3(b)). In the T2-chaotic 
regime it decays over the Ti-intervals, but it is magnified over the T2-intervals. The behavior 
is shown in Fig. E(c). 

Synchronization. The sub-LE Ao and the Ti-spectrum Ai appear naturally in the context 
of synchronization in a network of networks. In this case the interaction delays T\ within 
a subnetwork are much shorter than the connection delays r 2 between the different subnet- 
works. We illustrate this with a simple hierarchical network of dynamical units, described 
by 

i 

+e(l-^ B jT ) Cm S f(x k t LJ- (9) 

sk 

The coupling topology within the m-th subnetwork is described by the matrix A^ m \ the 
matrix C describes the coupling architecture of the networks, and the matrix B^ 1 ^ models 
the coupling between the m-th and the s-th subnetwork. We assume that all the elements 
receive the same amount of input, both from within and from outside the subnetwork (i.e. 
all the matrices have a row sum equal to 1), so that all the nodes can synchronize identi- 
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cally. Moreover, we consider identical rows for the £>-matrices, so that all the elements of a 
subnetwork receive the same input from outside. 

To determine the stability of a (cluster) synchronized state s™, the network model Eq. 
Q is linearized along the corresponding synchronization manifold. Evaluating this linear 
system along the transverse eigendirections yields 

5x t+1 = (1 - e)f{s?)5x t + ai m W'(*r ri )^_ T1 , (10) 

with a^ 1 the transverse eigenvalues of the connection matrix A^ m \ Integration of Eq. ( fit)]) 



reveals the master stability function \{ga) [27|. The explicit dependence on the long delay 
connections vanishes, as all the nodes receive the same external input. Hence, Eq.f llOp and 
the corres pon ding master stability function take the same form as for a network with a single 



delay |19l . |28j. The difference lies within the dynamics of the synchronized state s™, which 
now also depends on the elements outside the subnetwork and the connection delay t^- In 
the strongly chaotic regime, the largest transverse LE is approximately A > 0, so that the 
elements cannot synchronize. In the Ti-chaotic regime, identical or cluster synchronization 
in the m-th subnetwork is stable if \<J^\ < e _AmaxT1 



19j; the synchronization pattern in 



the subnetwork thus depends on the coupling topology A^ m \ For T2-chaotic behavior all 
the nodes of a subnetwork synchronize completely irrespective of the coupling architecture: 
the elements of the subnetwork show consistent behavior with respect to the common input 
from the rest of the network. 

The stability of identical synchronization of the full network is then governed by the 
following equation 

Sx t+ i = (1- e)f'(st)Sxt + eKf'(st- Tl )Sxt- T1 

+a c e{l - K)f(s t _ T2 )5x t _ T2 , (11) 

with oc the transverse eigenvalues of the intra-network coupling matrix C. Synchronization 
between subnetworks is only possible in the T2-chaotic regime in the limit of long delays, on 
the condition that \a c \ < e~ AmaxT2 holds. 

As an example we consider a network of four globally coupled subnetworks, coupled 
through their mean fields (Fig. HJ^a)). The subnetworks are bidirectional rings of four 
elements. As individual dynamics we choose again tent maps (Eq. flB])). We increase 
the total coupling strength e, such that the system undergoes a transition from strong 
to 7~i - to r 2 -chaos. Figure H] shows the crosscorrelations between several different network 
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FIG. 4: (Color online) Panel (a) shows a sketch of the network topology. In panel (b) the sub- 
LE Ao and Ai jmax and the maximal LE A max are plotted for a hierarchical network of tent maps 
(described by Eqs. Q and d5|)) as a function of the total coupling strength e. The dashed lines 
indicate the synchronization transitions. Panel (c) shows the crosscorrelations between nodes A 
and B (red dots), B and C (green squares) and A and D (blue crosses) for varying total coupling 
strength e. Parameters are a = 0.4, k = 0.5, t\ = 50, and t<i = 500. 

elements as function of e, together with the different sub-LEs Ao and Ai jmax . The sublattice 
synchronization between the diagonal elements B and C in a subnetwork is governed by 
Eq. ( |T0l) . with a transverse eigenvalue a a = 0. The nodes B and C thus synchronize when 
Ao < 0. The small difference between the two transition points is caused by numerical 
inaccuracy. The synchronization between all elements in a subnetwork is governed by the 
same equation (Eq. (llOp ). but the transverse eigenvalue with maximal magnitude is given 
by oa = —1- Consequently, the nodes A and B synchronize when Ai < 0. The inset in Fig. 
S] shows an exact agreement between these two points. We find that the whole network, and 
thus the nodes A and D, synchronizes for a slightly higher coupling strength. 

In conclusion, we showed that a hierarchy of time scales emerges in systems with several 
delays. These time scales can be characterized by the different components in the spectrum 
of LEs. Depending on the leading components of the spectrum, one can distinguish between 
strong chaos or r^-chaos. In the r^-chaotic regime, small perturbations evolve on the time 
scale of the time delay r^. Although these results are relevant for any systems with different 
and well-separated time delays, especially interesting is the application to the network of 
networks, where time delays within a subnetwork are shorter than the corresponding time 
delays between the different subnetworks. We showed, that in such a case, the units within 
a subnetwork can only synchronize when strong chaos is absent and the maximal LE scales 
either with the shorter or with the longer delay. The total synchronization of all elements 
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involving also all subnetworks is, however, possible only when the whole LE spectrum scales 
with the longer delay. 
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